% setup env
clear all;
close all;
clc;

% schem parameter
R = 120;
L = 100e-6;
C = 33e-12;

% Resonance friquecy
f0 = 1 / (2*pi * sqrt(L*C));
fprintf('f0 = %f MHz\n', f0/1e6); 
%Discrete time
T = 10e-5 / f0;
Tmod = 4 / f0;

% time vector
t = 0:T:Tmod;

% init array
U = nan(1, length(t));
dU = nan(1, length(t));
I = nan(1, length(t));

% initial condition
U(1) = 0;
dU(1) = 0;
I(1) = 0;

%Evaluate AFC
Fmin = 2/Tmod;
Fmax = 2 * f0;
A = 1;
% Friquency array
f = Fmin:((Fmax-Fmin) / 1e+2):Fmax;
K = nan(1, length(f));

h = waitbar(0);
for i = 1:length(f)
    E = A * cos(2*pi*f(i)*t);
    
    for k = 2:length(t)
        U(k) = U(k-1) + dU(k-1)*T;
        I(k) = I(k-1) + T*(E(k) - I(k-1)*R -U(k-1)) / L;
        dU(k) = I(k) / C;
    end
    Us = U(fix(end/2):end);
    K(i) = (max(Us)-min(Us)) / 2;
    waitbar(i/length(f),h);
end

% plottig
graphics_toolkit gnuplot
figure;
[y1, x1] = max(K);
plot(f,K,'linewidth',2,[f(x1)], [y1], '*');
grid on
text (f(x1) , y1, sprintf("x = %e\ny = %e" ,f(x1), y1))
legend({'K(f)'}, 'location', 'northeast');
title('Amplitude-Friquency Characteristic');
xlabel('f, Hz');
ylabel('K');

print -deps -color "./AFC.eps"
print -dpng -color "./AFC.png"
